# Perform all necessary regressions
theData = get(dataset)

outList = list()
dvtype = c("pred", "count")
themes = c("INF", "ANIM", "BELLIG", "IRR")

coefList = list()
for (i in 1:length(dvtype)) {
  
  themeFits = list()
  for (t in 1:length(themes)) {
    
    # If probability measure, use quasibinomial model
    if (dvtype[i]=="pred") {
      thedv = paste0(dvtype[i], themes[t])
      
      oneIVfits = list()
      
      fit2 = glm(paste(thedv, theIV, predControls, sep=" "), data=theData, family=quasibinomial,
                 control = glm.control(epsilon=0.00001, maxit=100, trace=FALSE))
      fit2c = coeftest(fit2, vcov=vcovCL, cluster=~ccode)
      
      # If count measure, use Poisson model
    } else if (dvtype[i]=="count") {
      
      thedv = paste0(dvtype[i], themes[t])
      
      fit2 = glm(paste(thedv, theIV, countControls, sep=" "), data=theData, family="poisson",
                 control = glm.control(epsilon=0.00001, maxit=100, trace=FALSE))
      fit2c = coeftest(fit2, vcov=vcovCL, cluster=~ccode)
      
    }
    
    # Store individual results
    oneIVfits[[1]] = fit2
    oneIVfits[[2]] = fit2c
    
    # Label the results
    names(oneIVfits) = c("gs_full", "gs_full_cluster")    
    
    # Store all results
    themeFits[[t]] = oneIVfits
    
  }
  names(themeFits) = themes
  
  coefList[[i]] = themeFits
}
names(coefList) = dvtype


## Results for the probabilities 
p_inf_full = coefList$pred$INF$gs_full
p_anim_full = coefList$pred$ANIM$gs_full
p_bel_full = coefList$pred$BELLIG$gs_full
p_irr_full = coefList$pred$IRR$gs_full

p_inf_full_c = coefList$pred$INF$gs_full_cluster
p_anim_full_c = coefList$pred$ANIM$gs_full_cluster
p_bel_full_c = coefList$pred$BELLIG$gs_full_cluster
p_irr_full_c = coefList$pred$IRR$gs_full_cluster

## Results for the counts 
c_inf_full = coefList$count$INF$gs_full
c_anim_full = coefList$count$ANIM$gs_full
c_bel_full = coefList$count$BELLIG$gs_full
c_irr_full = coefList$count$IRR$gs_full

c_inf_full_c = coefList$count$INF$gs_full_cluster
c_anim_full_c = coefList$count$ANIM$gs_full_cluster
c_bel_full_c = coefList$count$BELLIG$gs_full_cluster
c_irr_full_c = coefList$count$IRR$gs_full_cluster

# Save the models themselves
outList[[1]] = p_inf_full
outList[[2]] = p_anim_full
outList[[3]] = p_bel_full
outList[[4]] = p_irr_full
outList[[5]] = c_inf_full
outList[[6]] = c_anim_full
outList[[7]] = c_bel_full
outList[[8]] = c_irr_full

outList[[9]] = p_inf_full_c
outList[[10]] = p_anim_full_c
outList[[11]] = p_bel_full_c
outList[[12]] = p_irr_full_c
outList[[13]] = c_inf_full_c
outList[[14]] = c_anim_full_c
outList[[15]] = c_bel_full_c
outList[[16]] = c_irr_full_c

# Add labels to items in the list
names(outList) = c("p_inf", "p_anim", "p_bel", "p_irr",
                   "c_inf", "c_anim", "c_bel", "c_irr",
                   "p_inf_c", "p_anim_c", "p_bel_c", "p_irr_c",
                   "c_inf_c", "c_anim_c", "c_bel_c", "c_irr_c")


#### Make a prediction function ####
makePredict = function(theModel, theVar, datasetP) {
  
  # The dataset
  theData = datasetP
  
  # Create a mean observation
  meanob = theData %>% dplyr::select(predINF, predANIM, predBELLIG, predIRR, 
                                     countINF, countANIM, countBELLIG, countIRR,
                                     logYSI, conf_high, personal, leaderTenure, UStrade, USmilaid, logWords, #nWords,
                                     globalsouth, democracy, leaderMention, USdefense, contains("Topic"))
  meanob = data.frame(t(apply(meanob, 2, mean, na.rm=T)))
  meanob$admin = "Johnson"
  meanob$admin = factor(meanob$admin, levels=c("Kennedy", "Johnson", "Nixon", "Ford"))
  meanob$date = mean(theData$date)
  
  # Duplicate the mean observation
  meanobs = rbind(meanob, meanob)
  meanobs$country = "India"
  
  # Change "theVar" to assume minimum and maximum values
  meanobs[,theVar] = c(quantile(theData[,..theVar], 0.00, na.rm=T), quantile(theData[,..theVar], 1, na.rm=T))
  
  # Use these two observations to create predictions
  thePred = predict(theModel, newdata=meanobs, type="response")
  return(thePred)
}


# Make predictions using the function above
createPredictions = function(outputs) {
  
  p_inf_full = outputs[[1]]
  p_anim_full = outputs[[2]]
  p_bel_full = outputs[[3]]
  p_irr_full = outputs[[4]]
  c_inf_full = outputs[[5]]
  c_anim_full = outputs[[6]]
  c_bel_full = outputs[[7]]
  c_irr_full = outputs[[8]]
  
  dataset = p_inf_full[["data"]]
  
  suppressWarnings(
    if (grepl("globalsouth", theIV)==TRUE) {
      # Make predicted values
      p_inf_gs = makePredict(p_inf_full, "globalsouth", datasetP = dataset)
      p_inf_ysi = makePredict(p_inf_full, "logYSI", datasetP = dataset)
      p_anim_gs = makePredict(p_anim_full, "globalsouth", datasetP = dataset)
      p_anim_ysi = makePredict(p_anim_full, "logYSI", datasetP = dataset)
      p_bel_gs = makePredict(p_bel_full, "globalsouth", datasetP = dataset)
      p_bel_ysi = makePredict(p_bel_full, "logYSI", datasetP = dataset)
      p_irr_gs = makePredict(p_irr_full, "globalsouth", datasetP = dataset)
      p_irr_ysi = makePredict(p_irr_full, "logYSI", datasetP = dataset)
      
      c_inf_gs = makePredict(c_inf_full, "globalsouth", datasetP = dataset)
      c_inf_ysi = makePredict(c_inf_full, "logYSI", datasetP = dataset)
      c_anim_gs = makePredict(c_anim_full, "globalsouth", datasetP = dataset)
      c_anim_ysi = makePredict(c_anim_full, "logYSI", datasetP = dataset)
      c_bel_gs = makePredict(c_bel_full, "globalsouth", datasetP = dataset)
      c_bel_ysi = makePredict(c_bel_full, "logYSI", datasetP = dataset)
      c_irr_gs = makePredict(c_irr_full, "globalsouth", datasetP = dataset)
      c_irr_ysi = makePredict(c_irr_full, "logYSI", datasetP = dataset)
      
      # Put together predicted values
      p_inf = c(p_inf_gs, p_inf_ysi)
      p_anim = c(p_anim_gs, p_anim_ysi)
      p_bel = c(p_bel_gs, p_bel_ysi)
      p_irr = c(p_irr_gs, p_irr_ysi)
      
      c_inf = c(c_inf_gs, c_inf_ysi)
      c_anim = c(c_anim_gs, c_anim_ysi)
      c_bel = c(c_bel_gs, c_bel_ysi)
      c_irr = c(c_irr_gs, c_irr_ysi)
      
      predTable = rbind(p_inf, p_anim,
                        p_bel, p_irr, 
                        c_inf, c_anim,
                        c_bel, c_irr)
      print(round(predTable, 3))
      
    }
  )
}
